clear;
clc;

xa = -1.0;
xb = 1.0;
ya = -1.0;
yb = 1.0;
za = -1.0;
zb = 1.0;

vector_N = [2 4 8 16 32];

basis_type = "P1";
N_gauss_int2d = 9;
boundary_type = "Dirichlet";

Nx = vector_N(1);
Ny = Nx;
Nz = Nx;

N = reference_basis_function_3d(basis_type);

% plot_basis_function_3d(basis_type, 1);